{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_This setup code is required to run in an IPython notebook_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn\n",
    "\n",
    "seaborn.set_style('darkgrid')\n",
    "plt.rc(\"figure\", figsize=(16, 6))\n",
    "plt.rc(\"savefig\", dpi=90)\n",
    "plt.rc(\"font\",family=\"sans-serif\")\n",
    "plt.rc(\"font\",size=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sharpe Ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Sharpe Ratio is an important measure of return per unit of risk.  The example shows how to estimate the variance of the Sharpe Ratio and how to construct confidence intervals for the Sharpe Ratio using a long series of U.S. equity data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import arch.data.frenchdata\n",
    "\n",
    "ff = arch.data.frenchdata.load()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data set contains the Fama-French factors, including the excess market return."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            Mkt-RF          SMB          HML           RF\n",
      "count  1109.000000  1109.000000  1109.000000  1109.000000\n",
      "mean      0.659946     0.206555     0.368864     0.274220\n",
      "std       5.327524     3.191132     3.482352     0.253377\n",
      "min     -29.130000   -16.870000   -13.280000    -0.060000\n",
      "25%      -1.970000    -1.560000    -1.320000     0.030000\n",
      "50%       1.020000     0.070000     0.140000     0.230000\n",
      "75%       3.610000     1.730000     1.740000     0.430000\n",
      "max      38.850000    36.700000    35.460000     1.350000\n"
     ]
    }
   ],
   "source": [
    "excess_market = ff.iloc[:, 0]\n",
    "print(ff.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to construct a function that computes the Sharpe Ratio.  This function also return the annualized mean and annualized standard deviation which will allow the covariance matrix of these parameters to be estimated using the bootstrap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sharpe_ratio(x):\n",
    "    mu, sigma = 12 * x.mean(), np.sqrt(12 * x.var())\n",
    "    values = np.array([mu, sigma, mu / sigma]).squeeze()\n",
    "    index = ['mu', 'sigma', 'SR']\n",
    "    return pd.Series(values, index=index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function can be called directly on the data to show full sample estimates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mu        7.919351\n",
       "sigma    18.455084\n",
       "SR        0.429115\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = sharpe_ratio(excess_market)\n",
    "params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproducibility\n",
    "\n",
    "All bootstraps accept the keyword argument `random_state` which can contain a NumPy `RandomState` instance. This allows the same pseudo random values to be used across multiple runs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Warning_\n",
    "\n",
    "_The bootstrap chosen must be appropriate for the data.  Squared returns are serially correlated, and so a time-series bootstrap is required._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bootstraps are initialized with any bootstrap specific parameters and the data to be used in the bootstrap.  Here the `12` is the average window length in the Stationary Bootstrap, and the next input is the data to be bootstrapped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6wAAAF9CAYAAAAAxuXnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5RWdb0/8Pf4cBeQ2yN5+4XokfIC4iVEJRJvJXYSTY8rl5qmZiddLTEQLbMUuaiZHVS8pFamlkuNUjxax+osL53jJTpL1DA9kYoZA4gKAwjD/P5oOacJ0AH2zOwZXq+1Zq159t7PZ3/2Wl/27Df7VtPQ0NAQAAAAKJmt2roBAAAAWB+BFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKKVObd0AAGxJ5s2bl+uvvz7PPvtsamtrs80222TQoEEZMWJEzj333CTJySefnCeffLLxO127ds2HP/zhHHfccTnllFOy1Vb+vxmALYPACgCt5Jlnnsmpp56abbfdNuPGjcuHPvSh/PWvf83//M//ZObMmY2BNUmq1WomTJiQJHnzzTfzwAMPZOrUqVmyZEnGjx/fVpsAAK2qpqGhoaGtmwCALcFZZ52V3//+93nooYfSr1+/JvP++te/ZuDAgUn+doa1trY2Dz30UOP8VatW5VOf+lSWLl2ap556KpVKpVV7B4C24JoiAGglr7zySnbdddd1wmqSxrC6IV27ds2ee+6Z5cuXZ/HixS3VIgCUisAKAK1khx12yAsvvJA//OEPm/T9BQsWpKamJr179y64MwAoJ/ewAkArOeOMM3L66adn3Lhx2XPPPbPffvtlxIgRGTlyZLp27dpk2bVr12bJkiVJkqVLl+aee+7J3Llzc+ihh6Zbt25t0T4AtDr3sAJAK/rd736X733ve/ntb3+burq6JEnPnj1z0UUX5bjjjkuy7lOC33PkkUfmsssuyzbbbNOqPQNAW3GGFQBa0T777JPrr78+9fX1mTdvXn7961/ntttuy0UXXZTtt98+I0eOTJJst912mTJlStauXZs///nPufHGG7Nw4cJ07ty5jbcAAFqPe1gBoA1UKpXsvvvu+fKXv5wZM2YkSX7+8583zu/WrVsOPPDAHHzwwTnppJNy880359lnn83VV1/dVi0DQKsTWAGgjQ0dOjRJsnDhwg0uM2TIkBxzzDH58Y9/nAULFrRWawDQpgRWAGglv/3tb7N27dp1pv/nf/5nkmTw4MHv+/0vfOELWbNmTW655ZYW6Q8AysY9rADQSi6//PLU1dXlsMMOyy677JK1a9fm+eefz89+9rP06dMnp5566vt+f/DgwfnEJz6Re++9N1/+8pfTv3//VuocANqGM6wA0EomTpyYAw88MI899limT5+eyy+/PI8//ng+/elP55577smOO+74gTW+8IUvZOXKlfn+97/f8g0DQBvzWhsAAABKyRlWAAAASklgBQAAoJQEVgAAAEpJYAUAAKCUBFYAAABKqV28h3Xt2rWpr9/8hxlXKjWF1IG2YPzSXhm7tGfGL+2Z8Ut70rlzZb3T20Vgra9vyNKldZtdp0+fHoXUgbZg/NJeGbu0Z8Yv7ZnxS3tSrfZa73SXBAMAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKXUrMB644035rjjjss+++yTAw44IGeffXZefPHFJss0NDRkxowZOfjggzN06NCcfPLJ+eMf/9hkmbfeeisTJkzIvvvum3333TcTJkzI22+/XdzWAAAA0GE0K7A++eST+dznPpcf//jH+cEPfpBKpZLTTjstS5cubVzm5ptvzq233pqLL74499xzT/r165fTTjsty5Yta1zm/PPPz/PPP5+bb7453/ve9/L8889n4sSJxW8VAAAA7V5NQ0NDw8Z+afny5dlvv/1y3XXXZcyYMWloaMioUaNy0kkn5Utf+lKSZOXKlRk5cmQuuOCCnHjiiXn55Zdz1FFH5c4778y+++6bJHn66adz0kkn5d///d8zePDgDa5v9er6Ql567OXJtGfGL+2VsUt7ZvzSnhm/tCfVaq/1Tt+ke1iXL1+etWvXpnfv3kmS1157LbW1tTnooIMal+nWrVv233//zJkzJ0kyZ86c9OjRI/vss0/jMvvuu2969OjRuAwAAAC8p9OmfOnyyy/PRz/60QwfPjxJUltbmyQZMGBAk+X69++fhQsXJkkWLVqUfv36paampnF+TU1N+vXrl0WLFr3v+iqVmvTp02NTWv2HOlsVUgfagvFLe2Xs0p4Zv7Rnxi8dwUYH1qlTp+aZZ57JXXfdlUql0mTe34fR9Vnf/IaGhg/8Xn19g0uC2eIZv7RXxi7tmfFLe2b80p4UcknwlClTMnv27PzgBz/ITjvt9HfFq0n+70zrexYvXtx41nXAgAFZvHhx/v6W2YaGhrz55pvp37//xrQBAADAFqDZZ1gnT56cBx98MLfffnt22WWXJvN23HHHVKvVPPHEExk6dGiSZNWqVXn66acbnwI8fPjw1NXVZc6cOY33sc6ZMyd1dXWNlxYDwD/q2bt7unfdpDtYNmjFqjVZ9vaKQmsCAMVr1hHAt771rfzsZz/Lddddl969ezeeSe3Ro0e23nrr1NTU5JRTTskNN9yQwYMHZ9CgQZk5c2Z69OiRo48+Okmyyy67ZNSoUbnkkkty2WWXpaGhIZdcckkOOeSQ931CMABbtu5dO2XQpNmF1pw/bWyWffBiAEAba1ZgvfPOO5Mkn//855tMP+ecc3LuuecmSc4888ysWrUql156ad56660MGzYst956a3r27Nm4/FVXXZXJkyfn9NNPT5KMGTMm3/jGN4rYDgAAADqYTXoPa2vzHlYwfmm/NnfsVqu9WuQMa23tO4XWpGOy76U9M35pTwp9DysAAAC0NIEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBS6tTWDQDQtnr27p7uXYv7c7Bi1Zose3tFYfVawsrV9alWexVWrz1sMwC0RwIrwBaue9dOGTRpdmH15k8bm2WFVWsZ3TpXtrhtBoD2yCXBAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCl1as5CTz31VG655ZY899xzWbhwYaZOnZpjjz22cf6QIUPW+73Pfe5zueSSS5IkkyZNyk9/+tMm84cNG5a77757U3sHAACgA2tWYK2rq8tuu+2WY445JhdccME68x977LEmn+fOnZuzzz47n/rUp5pMP/DAA3PFFVc0fu7cufOm9AwAAMAWoFmBdfTo0Rk9enSS5MILL1xnfrVabfL5kUceyaBBg/Kxj32syfQuXbqssywAAACsT+H3sC5fvjyzZ8/OCSecsM68Z555JiNHjsyRRx6Zr3/961m8eHHRqwcAAKCDaNYZ1o3xwAMPZPXq1Rk3blyT6aNGjcrhhx+eHXfcMQsWLMg111yTU089Nffdd1+6dOnyvjUrlZr06dNjs3urVLYqpA60BeOX9uTvx+qWMHZXrq5Ptdqr0HqVwqqxObaE8UvHZfzSERQeWO++++4ceuih6devX5PpY8eObfx9yJAh2WOPPTJmzJj85je/yRFHHPG+NevrG7J0ad1m99anT49C6kBbMH5pKUUGrff8/Vjd3LHbEv0VrVvnSgZNml1YvfnTxqa29p3C6rHp7Htpz4xf2pMN/b0v9JLgF154IXPnzl3v5cD/aODAgRk4cGDmz59fZAsAAAB0EIUG1p/85CfZYYcdcuCBB37gskuWLMnChQuz7bbbFtkCAAAAHUSzLglevnx5XnnllSTJ2rVr8/rrr+eFF17INttsk+233z5JsmLFitx///0544wzUlNTs873r7322hxxxBGpVqtZsGBBrr766vTr1y+HHXZYwZsEAABAR9CswDp37tyccsopjZ9nzJiRGTNmZNy4cZk2bVqS5MEHH8yKFSty7LHHrvP9SqWSF198MbNmzco777yTarWaESNG5JprrknPnj0L2hQAymB9DyBqD/ehAgDl06zAOmLEiMybN+99lznuuONy3HHHrXdet27dcsstt2x8dwC0Oy3xACIAYMtU+HtYAQAAoAgCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApdSprRsAYOP07N093bvafQMAHZ8jHoB2pnvXThk0aXZh9eZPG1tYLQCAIrkkGAAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKKVObd0AANDUytX1qVZ7FVpzxao1Wfb2ikJrAkBLE1gBoGS6da5k0KTZhdacP21slhVaEQBankuCAQAAKKVmBdannnoqZ599dkaNGpUhQ4bkvvvuazJ/0qRJGTJkSJOfE044ocky7777bi677LKMGDEie++9d84+++y88cYbxW0JAAAAHUqzAmtdXV122223fO1rX0u3bt3Wu8yBBx6Yxx57rPHnpptuajL/8ssvz8MPP5yrr746d9xxR5YvX54vfvGLqa+v3/ytAAAAoMNp1j2so0ePzujRo5MkF1544XqX6dKlS6rV6nrnvfPOO7n33nszZcqUHHTQQUmSK664IoccckieeOKJjBo1alN6BwAAoAMr7KFLzzzzTEaOHJnevXtn//33z3nnnZf+/fsnSebOnZvVq1fn4IMPblx+u+22yy677JI5c+Z8YGCtVGrSp0+Pze6xUtmqkDrQFoxfYHPZh2w8+17aM+OXjqCQwDpq1Kgcfvjh2XHHHbNgwYJcc801OfXUU3PfffelS5cuWbRoUSqVSvr27dvke/3798+iRYs+sH59fUOWLq3b7D779OlRSB1oC8Yv7yn6dSdsOexDNp59L+2Z8Ut7sqHjm0IC69ixYxt/HzJkSPbYY4+MGTMmv/nNb3LEEUds8HsNDQ1FrB4AAIAOqEVeazNw4MAMHDgw8+fPT5IMGDAg9fX1efPNN5sst2TJkgwYMKAlWgAAAKCda5HAumTJkixcuDDbbrttkmTPPfdM586d8/jjjzcu88Ybb+Tll1/O8OHDW6IFAAAA2rlmXRK8fPnyvPLKK0mStWvX5vXXX88LL7yQbbbZJttss02uvfbaHHHEEalWq1mwYEGuvvrq9OvXL4cddliSpFevXjnuuONyxRVXpH///unTp0+mTp2aIUOG5MADD2y5rQMAAKDdalZgnTt3bk455ZTGzzNmzMiMGTMybty4fPOb38yLL76YWbNm5Z133km1Ws2IESNyzTXXpGfPno3fueiii9KpU6ecd955WblyZUaOHJkrrrgilUql+K0CAACg3WtWYB0xYkTmzZu3wfm33HLLB9bo2rVrLr744lx88cXN7w4AKMTK1fWFPmF6xao1Wfb2isLqAcD6FPYeVgCgvLp1rmTQpNmF1Zs/bWyWFVYNANavRR66BAAAAJtLYAUAAKCUBFYAAABKSWAFAACglARWAAAASklgBQAAoJQEVgAAAEpJYAUAAKCUBFYAAABKqVNbNwDQkfXs3T3du9rVAgBsCkdRAC2oe9dOGTRpdqE1508bW2g9AICyckkwAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKzQqsTz31VM4+++yMGjUqQ4YMyX333dc4b/Xq1bnyyivz6U9/OnvvvXcOPvjgnH/++Xn99deb1Dj55JMzZMiQJj/nnXdesVsDAABAh9GpOQvV1dVlt912yzHHHJMLLrigybyVK1fm+eefz5e+9KV85CMfybJlyzJt2rScccYZ+fnPf55Onf5vFccee2zGjx/f+Llbt24FbQYAAAAdTbMC6+jRozN69OgkyYUXXthkXq9evXLbbbc1mXbppZdm7NixefnllzNkyJDG6d27d0+1Wt3cngEAANgCtMg9rMuWLUuSbLPNNk2mz549OyNGjMjYsWMzffr0xuUAAADgHzXrDOvGePfddzNt2rQccsgh+dCHPtQ4/eijj87222+fbbfdNi+99FK+/e1v5w9/+MM6Z2fXp1KpSZ8+PTa7t0plq0LqQFswfoGy2RL2Sfa9tGfGLx1BoYF1zZo1mTBhQt55553MnDmzybx/+Zd/afx9yJAh2WmnnXL88cfnueeeyx577PG+devrG7J0ad1m99enT49C6kBbMH7bp2q1V1u3AC1mS9gn2ffSnhm/tCcbOmYq7JLgNWvWZPz48Zk3b16+//3vp2/fvu+7/J577plKpZI///nPRbUAAABAB1LIGdbVq1dn/PjxefHFF3P77bc368FKL774Yurr6z2ECQAAgPVqVmBdvnx5XnnllSTJ2rVr8/rrr+eFF17INttsk2233TZf+cpX8uyzz+aGG25ITU1Namtrk/ztCcLdunXLK6+8kp///OcZPXp0+vbtm5dffjnTpk3L7rvvnn322afltg4AAIB2q1mBde7cuTnllFMaP8+YMSMzZszIuHHjcs455+SRRx5J8rf3rP69qVOn5thjj03nzp3zX//1X7n99tuzfPnybLfddhk9enTOOeecVCqVAjcHAACAjqJZgXXEiBGZN2/eBue/37wk2W677fKjH/1o4zoDAABgi9Yi72EFAACAzVX4e1gB2rOevbune1e7RgCAMnBUBvB3unftlEGTZhdWb/60sYXVAgDY0rgkGAAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBS6tTWDQAA7c/K1fWpVnsVVm/FqjVZ9vaKwuoB0DEIrADARuvWuZJBk2YXVm/+tLFZVlg1ADoKlwQDAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApdSswPrUU0/l7LPPzqhRozJkyJDcd999TeY3NDRkxowZOfjggzN06NCcfPLJ+eMf/9hkmbfeeisTJkzIvvvum3333TcTJkzI22+/XdyWAAAA0KF0as5CdXV12W233XLMMcfkggsuWGf+zTffnFtvvTXTpk3LzjvvnOuuuy6nnXZaHnroofTs2TNJcv755+cvf/lLbr755tTU1OTrX/96Jk6cmBtuuKHYLQK2KD17d0/3rs3alQEA0M406yhv9OjRGT16dJLkwgsvbDKvoaEhP/zhD3PWWWflyCOPTJJMnz49I0eOzAMPPJATTzwxL7/8ch599NHceeed2WeffZIk3/rWt3LSSSflf//3fzN48OAitwnYgnTv2imDJs0urN78aWMLqwUAwObZ7HtYX3vttdTW1uaggw5qnNatW7fsv//+mTNnTpJkzpw56dGjR2NYTZJ99903PXr0aFwGAAAA/t5mX0dXW1ubJBkwYECT6f3798/ChQuTJIsWLUq/fv1SU1PTOL+mpib9+vXLokWLPnAdlUpN+vTpsbmtplLZqpA60BaMX6CjK+M+zr6X9sz4pSMo7Mavvw+jzZ3f0NDwgd9Lkvr6hixdWrfJvb2nT58ehdSBtmD8rl+12qutWwAKUvQ+rqh73LfaqpIkWbFqTZa9vWKz60FrcexAe7KhY7rN3otXq9UkfzvTut122zVOX7x4ceNZ1wEDBmTx4sVNAmpDQ0PefPPN9O/ff3NbAABYR0vc476ssGoANMdmB9Ydd9wx1Wo1TzzxRIYOHZokWbVqVZ5++ulMnDgxSTJ8+PDU1dVlzpw5jfexzpkzJ3V1dRk+fPjmtgAAtHMrV9e7YgKAdTQrsC5fvjyvvPJKkmTt2rV5/fXX88ILL2SbbbbJ9ttvn1NOOSU33HBDBg8enEGDBmXmzJnp0aNHjj766CTJLrvsklGjRuWSSy7JZZddloaGhlxyySU55JBDPCEYAEi3zpVCz4YmnvoN0BE0K7DOnTs3p5xySuPnGTNmZMaMGRk3blymTZuWM888M6tWrcqll16at956K8OGDcutt97a+A7WJLnqqqsyefLknH766UmSMWPG5Bvf+EbBmwMAAEBH0azAOmLEiMybN2+D82tqanLuuefm3HPP3eAyffr0yVVXXbXxHQIAALBF2uz3sAIAAEBLEFgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKKVORRQZM2ZMFixYsM700aNH56abbsqMGTNy7bXXNpk3YMCAPP7440WsHgAAgA6okMB6zz33pL6+vvFzbW1tjj322HzqU59qnLbzzjvn9ttvb/xcqVSKWDUAAAAdVCGBtV+/fk0+33PPPenZs2c++clP/t+KOnVKtVotYnUAAABsAQoJrH+voaEh99xzT/75n/853bt3b5z+6quvZtSoUencuXOGDRuW8ePHZ6eddip69QAAAHQQhQfWxx9/PK+99lqOP/74xmlDhw7N1KlTM3jw4CxZsiQzZ87MiSeemAceeCB9+/b9wJqVSk369Omx2b1VKlsVUgfagvEL0Pbsh2lPHDvQERQeWO++++7stdde+ehHP9o4bfTo0U2WGTZsWA477LDMmjUrp5122gfWrK9vyNKldZvdW58+PQqpA23B+F2/arVXW7cAbEHsh2lPHDvQnmzomK7Q19osXrw4v/rVr3LCCSe873Jbb711dt1118yfP7/I1QMAANCBFBpY77vvvnTu3DlHHXXU+y63atWq/OlPf/IQJgAAADaosEuC33vY0tixY9OzZ88m86ZPn55DDjkk2223XZYsWZLrr78+dXV1GTduXFGrBwAAoIMpLLD+93//d+bPn58rr7xynXlvvPFGxo8fn6VLl6Zv377Ze++9c/fdd2eHHXYoavUAAAB0MIUF1gMOOCDz5s1b77zvfOc7Ra0GAACALUSh97ACAABAUQRWAAAASklgBQAAoJQKu4cV4IP07N093bva7QAA0DyOHIFW071rpwyaNLvQmvOnjS20HgAA5eGSYAAAAEpJYAUAAKCUBFYAAABKSWAFAACglARWAAAASklgBQAAoJQEVgAAAErJe1iBDerZu3u6d7WbAACgbTgSBTaoe9dOGTRpdmH15k8bW1gtAAA6PoEVAKAZVq6uT7Xaq7B6K1atybK3VxRWD6AjElgBAJqhW+dK4VedLCusGkDH5KFLAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKRUSWGfMmJEhQ4Y0+TnooIMa5zc0NGTGjBk5+OCDM3To0Jx88sn54x//WMSqAQAA6KA6FVVo5513zu233974uVKpNP5+880359Zbb820adOy884757rrrstpp52Whx56KD179iyqBQAAADqQwi4J7tSpU6rVauNPv379kvzt7OoPf/jDnHXWWTnyyCOz2267Zfr06Vm+fHkeeOCBolYPAABAB1NYYH311VczatSojBkzJuedd15effXVJMlrr72W2traJpcId+vWLfvvv3/mzJlT1OoBAADoYAq5JHjo0KGZOnVqBg8enCVLlmTmzJk58cQT88ADD6S2tjZJMmDAgCbf6d+/fxYuXNis+pVKTfr06bHZfVYqWxVSB9qC8QvQ8div05IcO9ARFBJYR48e3eTzsGHDcthhh2XWrFkZNmxYkqSmpmaT69fXN2Tp0rrN6jH52x+FIupAW2iL8Vut9mrV9QFsaRyX0JIc+9KebOi4s0Vea7P11ltn1113zfz581OtVpOk8UzrexYvXrzOWVcAAAB4T4sE1lWrVuVPf/pTqtVqdtxxx1Sr1TzxxBNN5j/99NMZPnx4S6weAACADqCQS4KnT5+eQw45JNttt12WLFmS66+/PnV1dRk3blxqampyyimn5IYbbsjgwYMzaNCgzJw5Mz169MjRRx9dxOoBAADogAoJrG+88UbGjx+fpUuXpm/fvtl7771z9913Z4cddkiSnHnmmVm1alUuvfTSvPXWWxk2bFhuvfVW72AFAABggwoJrN/5znfed35NTU3OPffcnHvuuUWsDgAAgC1Ai9zDCgAAAJurkDOsQDn07N093bv6Zw0AQMfgyBY6kO5dO2XQpNmF1Zs/bWxhtQAAYGO5JBgAAIBSElgBAAAoJYEVAACAUnIPKwBAG1i5uj7Vaq/C6q1YtSbL3l5RWD2AMhBYAQDaQLfOlcIflLessGoA5eCSYAAAAEpJYAUAAKCUBFYAAABKSWAFAACglDx0CQCgAyj6qcOJJw8DbU9gBQDoAIp+6nDiycNA23NJMAAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApdSprRsAAKCcVq6uT7Xaq7B6K1atybK3VxRWD+j4BFYAANarW+dKBk2aXVi9+dPGZllh1YAtgUuCAQAAKCWBFQAAgFIq5JLgG2+8Mb/4xS/ypz/9KV26dMnee++d8ePHZ7fddmtcZtKkSfnpT3/a5HvDhg3L3XffXUQLAAAAdDCFBNYnn3wyn/vc57LXXnuloaEh//Zv/5bTTjsts2fPTp8+fbZNC+IAAAwtSURBVBqXO/DAA3PFFVc0fu7cuXMRqwcAAKADKiSw3nLLLU0+X3HFFdlvv/3yu9/9LmPGjGmc3qVLl1Sr1SJWCQAAQAfXIk8JXr58edauXZvevXs3mf7MM89k5MiR6d27d/bff/+cd9556d+/f0u0AAAAQDvXIoH18ssvz0c/+tEMHz68cdqoUaNy+OGHZ8cdd8yCBQtyzTXX5NRTT819992XLl26vG+9SqUmffr02Oy+KpWtCqkDbcH4BaAj8Les9Th2oCMoPLBOnTo1zzzzTO66665UKpXG6WPHjm38fciQIdljjz0yZsyY/OY3v8kRRxzxvjXr6xuydGndZvfWp0+PQupAW2jO+C3y5e4A0BIci7Uex760Jxs6ji00sE6ZMiUPPvhgfvCDH2SnnXZ632UHDhyYgQMHZv78+UW2AAAAQAdRWGCdPHlyHnzwwdx+++3ZZZddPnD5JUuWZOHChdl2222LagEAAIAOpJDA+q1vfSs/+9nPct1116V3796pra1NkvTo0SNbb711li9fnmuvvTZHHHFEqtVqFixYkKuvvjr9+vXLYYcdVkQLAAAAdDCFBNY777wzSfL5z3++yfRzzjkn5557biqVSl588cXMmjUr77zzTqrVakaMGJFrrrkmPXv2LKIFAAAAOphCAuu8efPed363bt3WeVcrAAAAvJ8Wea0N0Dw9e3dP967N/2foKcAAAGxJBFZoQ927dsqgSbMLqzd/2tgPXggAANqJrdq6AQAAAFgfgRUAAIBSElgBAAAoJfew0mFt7AONPsiKVWuy7O0VhdUDAADen8BKh9USDzRaVlg1AADgg7gkGAAAgFISWAEAACglgRUAAIBScg8rNNPK1fWpVnu1dRsAALDFEFihmbp1rhT6EKfkbw9yAgAA1s8lwQAAAJSSwAoAAEApCawAAACUkntYAQBoFUU/wHDFqjVZ9vaKwuoB5SOwAgDQKop+gOEfLvtk4U/wF4KhXARWAADapZZ6gv+yQisCm8M9rAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCl5rQ2bpGfv7unetdjh471nAEBbW7m6vtB3uzq+gc0jsLJJunft5L1nAECHU/S7XR3fwOZxSTAAAACl5AwrpVH0JTgAAB3Nxt6W9UHHVi5ZpuxaPbDecccdueWWW1JbW5t/+qd/ykUXXZT99tuvtdughFriEhwAgI6k6NuyXLJM2bXqJcEPPvhgpkyZkrPPPjuzZs3K8OHDc+aZZ+b1119vzTYAAABoB1r1DOttt92WcePG5YQTTkiSXHzxxXn00Udz11135fzzz2/NVgpX9FNzV66uT7fOldLWAwDgg5X9lqeW6K/sx7Ht4TLoorNFe9jmDWm1wPruu+/mueeey+mnn95k+kEHHZQ5c+a0VhstpiUuzyh7PQAA3l/Zb3kqur+kfRzHlv0yaJd+/5+ahoaGhtZY0V//+td8/OMfz49+9KPsv//+jdOvvfba3H///Xn44Ydbow0AAADaiVZ/rU1NTU2zpgEAALBla7XA2rdv31QqldTW1jaZvnjx4gwYMKC12gAAAKCdaLXA2qVLl+yxxx554oknmkx/4oknMnz48NZqAwAAgHaiVZ8SfNppp2XixIkZOnRo9tlnn9x1111ZuHBhTjzxxNZsAwAAgHagVQPrUUcdlTfffDMzZ87MwoULs9tuu+Wmm27KDjvs0JptAAAA0A602lOCAQAAYGO0+lOCAQAAoDk6VGC94447MmbMmOy111459thj8/TTT7/v8k8++WSOPfbY7LXXXjn00ENz1113tVKnsK6NGb+/+MUvcvrpp+eAAw7I8OHDc/zxx+eRRx5pxW7h/2zsvvc9Tz/9dHbfffccffTRLdwhbNjGjt9333033/3udzNmzJjsueee+cQnPpEf/vCHrdQtNLWx4/f+++/PZz7zmQwbNiwHHXRQvvrVr67zBg8omw4TWB988MFMmTIlZ599dmbNmpXhw4fnzDPPzOuvv77e5V999dWcddZZGT58eGbNmpUvfvGLmTx5ch5++OFW7hw2fvw++eSTOeCAA3LTTTdl1qxZGT16dM4555xmBwUoysaO3fe89dZbueCCCzJy5MhW6hTWtSnj9/zzz8+jjz6ayy67LA899FC++93vZsiQIa3YNfzNxo7fZ555JhMnTsy4cePywAMP5LrrrsvLL7+cr371q63cOWycDnMP6/HHH58hQ4Zk8uTJjdOOOOKIHHnkkTn//PPXWf7KK6/ML3/5y/ziF79onPa1r30tL730Un7yk5+0Ss/wno0dv+vz2c9+Nvvtt18mTZrUUm3COjZ17J5zzjn5yEc+koaGhjz88MN54IEHWqNdaGJjx+9jjz2Wr3zlK/nlL3+Zfv36tWarsI6NHb+33HJLfvSjH+XXv/5147R77703kydPzpw5c1qlZ9gUHeIM67vvvpvnnnsuBx10UJPpBx100Ab/Af7+979fZ/mDDz44c+fOzerVq1usV/hHmzJ+12f58uXp3bt30e3BBm3q2L3jjjuyaNGifOlLX2rpFmGDNmX8/sd//Ef22muvfP/738/HP/7xHHHEEZk8eXKWL1/eGi1Do00Zv/vss09qa2vzq1/9Kg0NDVmyZEkefPDBfPzjH2+NlmGTdYjA+uabb6a+vj4DBgxoMr1///4bvC5/0aJF6d+/f5NpAwYMyJo1a/Lmm2+2WK/wjzZl/P6jO+64I2+88UY+85nPtESLsF6bMnbnzZuX6667LldeeWUqlUprtAnrtSnj99VXX80zzzyTP/zhD5kxY0YuvvjiPProo7nwwgtbo2VotCnjd/jw4fn2t7+dr371q9lzzz0zcuTINDQ0ZPr06a3RMmyyDhFY31NTU9OsaRua997V0e/3HWgpGzt+3/Pwww/niiuuyFVXXeWdxrSJ5o7dd999N+PHj8/EiROz0047tUZr8IE2Zt/b0NCQmpqafPvb386wYcMyatSoXHzxxXn44YezaNGilm4V1rEx4/ell17K5MmT86//+q+59957873vfS+1tbX5xje+0dJtwmbp1NYNFKFv376pVCrr/I/S4sWL1/mfp/cMGDBgnT8uixcvTqdOndKnT58W6xX+0aaM3/c8/PDDmThxYqZPn55DDz20JduEdWzs2F24cGFeeumlXHTRRbnooouSJGvXrk1DQ0N233333HTTTTn44INbpXfYlH1vtVrNwIED06tXr8Zpu+yyS5Lk9ddf/8B9NhRlU8bvjTfemKFDh+aMM85IknzkIx9J9+7dc9JJJ+W8887Ldttt1+J9w6boEGdYu3Tpkj322CNPPPFEk+lPPPFEhg8fvt7v7L333utdfs8990znzp1brFf4R5syfpO/PR1wwoQJmTp1aj75yU+2dJuwjo0duwMHDsz999+fWbNmNf6ceOKJ+fCHP9z4hEtoLZuy791nn32ycOHCJveszp8/P0lc4UKr2pTxu3LlynVuxXjvcwd5BisdVOWb3/zmN9u6iSL07NkzM2bMSLVaTbdu3XL99dfn6aefzpQpU9K7d+9MnDgxv/zlL3P44YcnSf7f//t/ufnmm7N48eLssMMOeeSRR3LDDTdk0qRJ2XXXXdt4a9jSbOz4nT17diZOnJgJEybkE5/4ROrq6lJXV5fVq1enW7dubbw1bEk2ZuxWKpX079+/yc+zzz6b+fPn59xzz02XLl3aenPYwmzsvnfnnXfOvffem7lz52bXXXfN/Pnzc+mll2b//ffPcccd18Zbw5ZmY8fvypUrc8stt6Rv377p06dPXnrppVx++eWpVquNZ12hjDrEJcFJctRRR+XNN9/MzJkzs3Dhwuy222656aabGv/H8y9/+UuT5XfaaafcdNNNmTp1au66665su+22+drXvpYjjzyyLdpnC7ex4/fHP/5x1qxZkylTpmTKlCmN0z/2sY/l9ttvb9Xe2bJt7NiFMtnY8bv11lvntttuy+TJk/PZz342vXv3zmGHHdbs149BkTZ2/B577LFZvnx57rjjjkyfPj29evXKiBEjMmHChLZoH5qtw7yHFQAAgI6lQ9zDCgAAQMcjsAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAACl9P8BCc1XPoqx+KcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from arch.bootstrap import StationaryBootstrap\n",
    "# Initialize with entropy from random.org\n",
    "entropy = [877788388, 418255226, 989657335, 69307515]\n",
    "rs = np.random.RandomState(entropy)\n",
    "\n",
    "bs = StationaryBootstrap(12, excess_market, random_state=rs)\n",
    "results = bs.apply(sharpe_ratio, 2500)\n",
    "SR = pd.DataFrame(results[:, -1:], columns=['SR'])\n",
    "fig = SR.hist(bins=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             mu     sigma        SR\n",
      "mu     3.835234 -0.700250  0.222532\n",
      "sigma -0.700250  3.517844 -0.119008\n",
      "SR     0.222532 -0.119008  0.014912\n",
      "\n",
      "\n",
      "mu       1.958375\n",
      "sigma    1.875592\n",
      "SR       0.122114\n",
      "Name: Std Errors, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "cov = bs.cov(sharpe_ratio, 1000)\n",
    "cov = pd.DataFrame(cov, index=params.index, columns=params.index)\n",
    "print(cov)\n",
    "se = pd.Series(np.sqrt(np.diag(cov)), index=params.index)\n",
    "se.name = 'Std Errors'\n",
    "print('\\n')\n",
    "print(se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              mu      sigma        SR\n",
      "Lower   4.319805  14.559092  0.186293\n",
      "Upper  12.004185  21.517835  0.660455\n"
     ]
    }
   ],
   "source": [
    "ci = bs.conf_int(sharpe_ratio, 1000, method='basic')\n",
    "ci = pd.DataFrame(ci, index=['Lower', 'Upper'], columns=params.index)\n",
    "print(ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative confidence intervals can be computed using a variety of methods.  Setting `reuse=True` allows the previous bootstrap results to be used when constructing confidence intervals using alternative methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              mu      sigma        SR\n",
      "Lower   3.834517  15.392333  0.197774\n",
      "Upper  11.518896  22.351076  0.671937\n"
     ]
    }
   ],
   "source": [
    "ci = bs.conf_int(sharpe_ratio, 1000, method='percentile', reuse=True)\n",
    "ci = pd.DataFrame(ci, index=['Lower', 'Upper'], columns=params.index)\n",
    "print(ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal Block Length Estimation\n",
    "\n",
    "The function `optimal_block_length` can be used to estimate the optimal block lengths for the Stationary and Circular bootstraps. Here we use the squared market return since the Sharpe ratio depends on the mean and the variance, and the autocorrelation in the squares is stronger than in the returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        stationary   circular\n",
      "Mkt-RF   47.766787  54.679322\n"
     ]
    }
   ],
   "source": [
    "from arch.bootstrap import optimal_block_length\n",
    "opt = optimal_block_length(excess_market ** 2)\n",
    "print(opt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can repeat the analysis above using the estimated optimal block length.  Here we see that the extremes appear to be slightly more extreme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6wAAAF9CAYAAAAAxuXnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5RWdb0/8Pc4XAYEHC6P5O2E6IFKBUENUYnEW4mdQtPjqiWmqdlJVksMRM0sRS5qZgcV71amlkuNUjxZp8tZJZ3jJTpL0jA9kYbHuInKDCAM8/ujn5NzuDjgnpnN8HqtNWvx7L2fz/7stb7sed6z9/PdVY2NjY0BAACAktmpvRsAAACATRFYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUurU3g0AwI5k4cKFufHGG/P0009n6dKl2WWXXTJgwICMGDEiEyZMSJKcdtppefzxx5ve07Vr17z3ve/NSSedlPHjx2ennfy9GYAdg8AKAG3kqaeeyumnn55dd90148aNy3ve85789a9/zX//939n9uzZTYE1SSqVSiZNmpQkefXVV/Pwww9n+vTpWbFiRSZOnNhehwAAbaqqsbGxsb2bAIAdwTnnnJPf/e53+fGPf5w+ffo0W/fXv/41/fv3T/K3K6xLly7Nj3/846b1a9euzUc/+tGsXLkyTzzxRKqrq9u0dwBoD+4pAoA28uKLL2bffffdKKwmaQqrm9O1a9fsv//+qaury/Lly1urRQAoFYEVANrIHnvskWeffTZ/+MMftun9ixcvTlVVVXr16lVwZwBQTr7DCgBt5KyzzsqZZ56ZcePGZf/998/BBx+cESNGZOTIkenatWuzbTds2JAVK1YkSVauXJn7778/CxYsyFFHHZWampr2aB8A2pzvsAJAG/rtb3+b2267Lb/5zW9SX1+fJOnRo0cuvvjinHTSSUk2niX4Lccdd1yuuOKK7LLLLm3aMwC0F1dYAaANDR8+PDfeeGMaGhqycOHC/OIXv8idd96Ziy++OLvvvntGjhyZJNltt90ybdq0bNiwIX/+859z8803Z8mSJencuXM7HwEAtB3fYQWAdlBdXZ0PfOAD+cIXvpBZs2YlSX70ox81ra+pqclhhx2WI444Ip/+9Kdz66235umnn861117bXi0DQJsTWAGgnQ0ZMiRJsmTJks1uM3jw4HziE5/I9773vSxevLitWgOAdiWwAkAb+c1vfpMNGzZstPw//uM/kiQDBw7c4vs/+9nPZv369bn99ttbpT8AKBvfYQWANnLllVemvr4+Rx99dPbZZ59s2LAhzzzzTH74wx+mtrY2p59++hbfP3DgwHz4wx/OAw88kC984Qvp27dvG3UOAO3DFVYAaCOTJ0/OYYcdll//+teZOXNmrrzyyjz22GP52Mc+lvvvvz977rnnO9b47Gc/mzVr1uRb3/pW6zcMAO3MY20AAAAoJVdYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAope3iOawbNmxIQ0NxkxlXV1cVWg/ai7FMR2Ac01EYy3QUxjLtoXPn6k0u3y4Ca0NDY1aurC+sXm1t90LrQXsxlukIjGM6CmOZjsJYpj1UKj03udwtwQAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUUosC680335yTTjopw4cPz6GHHppzzz03zz33XLNtGhsbM2vWrBxxxBEZMmRITjvttPzxj39sts1rr72WSZMm5aCDDspBBx2USZMm5fXXXy/uaAAAAOgwWhRYH3/88XzqU5/K9773vXz7299OdXV1zjjjjKxcubJpm1tvvTV33HFHLr300tx///3p06dPzjjjjKxatappmwsuuCDPPPNMbr311tx222155plnMnny5OKPCgAAgO1eVWNjY+PWvqmuri4HH3xwbrjhhowZMyaNjY0ZNWpUPv3pT+fzn/98kmTNmjUZOXJkLrzwwpx66ql54YUXcvzxx+eee+7JQQcdlCR58skn8+lPfzr/9m//loEDB252f+vWNWTlyvptPMSN1dZ2L7QetBdjmY7AOKajMJbpKIxl2kOl0nOTy7fpO6x1dXXZsGFDevXqlST5y1/+kqVLl+bwww9v2qampiaHHHJI5s+fnySZP39+unfvnuHDhzdtc9BBB6V79+5N2wAAAMBbOm3Lm6688sq8//3vz7Bhw5IkS5cuTZL069ev2XZ9+/bNkiVLkiTLli1Lnz59UlVV1bS+qqoqffr0ybJly7a4v+rqqtTWdt+WVjdTb6dC60F7MZbpCNpjHDckqelcXVi9NesaUlw1tlfOyXQUxjJlstWBdfr06Xnqqady7733prq6+a/nt4fRTdnU+sbGxnd8X0NDo1uCYROMZTqC9hjHlUrPDJgyt7B6i2aMzdKlbxRWj+2TczIdhbFMeyjkluBp06Zl7ty5+fa3v5299trrbcUrSf5+pfUty5cvb7rq2q9fvyxfvjxv/8psY2NjXn311fTt23dr2gAAAGAH0OIrrFOnTs0jjzySu+66K/vss0+zdXvuuWcqlUrmzZuXIUOGJEnWrl2bJ598smkW4GHDhqW+vj7z589v+h7r/PnzU19f33RrMQBtr0evbunWdZu+IbJJq9euz6rXVxdWDwDYcbXoE8rXvva1/PCHP8wNN9yQXr16NV1J7d69e3beeedUVVVl/PjxuemmmzJw4MAMGDAgs2fPTvfu3XPCCSckSfbZZ5+MGjUql112Wa644oo0Njbmsssuy5FHHrnFGYIBaF3dunYq/PbYVe+8GQDAO2pRYL3nnnuSJJ/5zGeaLT/vvPMyYcKEJMnZZ5+dtWvX5vLLL89rr72WoUOH5o477kiPHj2atr/mmmsyderUnHnmmUmSMWPG5Ctf+UoRxwEAAEAH06LAunDhwnfcpqqqKhMmTGgKsJtSW1uba665puXdAQAAsMPapuewAgAAQGsTWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAopU4t2eiJJ57I7bffnt///vdZsmRJpk+fnhNPPLFp/eDBgzf5vk996lO57LLLkiRTpkzJD37wg2brhw4dmvvuu29bewcAAKADa1Fgra+vz6BBg/KJT3wiF1544Ubrf/3rXzd7vWDBgpx77rn56Ec/2mz5YYcdlquuuqrpdefOnbelZwAAAHYALQqso0ePzujRo5MkF1100UbrK5VKs9c/+9nPMmDAgHzwgx9strxLly4bbQsAAACbUvh3WOvq6jJ37tyccsopG6176qmnMnLkyBx33HH58pe/nOXLlxe9ewAAADqIFl1h3RoPP/xw1q1bl3HjxjVbPmrUqBxzzDHZc889s3jx4lx33XU5/fTT8+CDD6ZLly5brFldXZXa2u6F9VhdvVOh9aC9GMuU1daMy44yjjvCMfDudJSxDMYyZVJ4YL3vvvty1FFHpU+fPs2Wjx07tunfgwcPzn777ZcxY8bkl7/8ZY499tgt1mxoaMzKlfWF9Vhb273QetBejGWKUKn0LLzm1ozL9hjH7X3MdEzOyXQUxjLtYXO/mwsNrM8++2wWLFiQiRMnvuO2/fv3T//+/bNo0aIiWwCgg+nRq1u6dS3876sAwHag0E8A3//+97PHHnvksMMOe8dtV6xYkSVLlmTXXXctsgWADm9HC3DdunbKgClzC625aMbYd94IAGh3LfrEU1dXlxdffDFJsmHDhrz88st59tlns8suu2T33XdPkqxevToPPfRQzjrrrFRVVW30/uuvvz7HHntsKpVKFi9enGuvvTZ9+vTJ0UcfXfAhAXRsRQc44Q0AKKsWBdYFCxZk/PjxTa9nzZqVWbNmZdy4cZkxY0aS5JFHHsnq1atz4oknbvT+6urqPPfcc5kzZ07eeOONVCqVjBgxItddd1169OhR0KEAAADQkbQosI4YMSILFy7c4jYnnXRSTjrppE2uq6mpye2337713QEAALDDKvw5rAAAAFAEgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKKUd58nzALSJNesaUqn03Kr3bO32AMCOQWAFoFA1naszYMrcwuotmjG2sFoAwPbFLcEAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSWYIB4F3alkf5vFO9ms7VhdVLktVr12fV66sLrQkArU1gBYB3qTUe5VNkvbdqriq0IgC0PrcEAwAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKUksAIAAFBKAisAAACl1Km9GwAAWt+adQ2pVHoWVm/12vVZ9frqwuoBwKYIrACwA6jpXJ0BU+YWVm/RjLFZVVg1ANg0twQDAABQSgIrAAAApdSiwPrEE0/k3HPPzahRozJ48OA8+OCDzdZPmTIlgwcPbvZzyimnNNvmzTffzBVXXJERI0bkwAMPzLnnnptXXnmluCMBAACgQ2lRYK2vr8+gQYNyySWXpKamZpPbHHbYYfn1r3/d9HPLLbc0W3/llVfm0UcfzbXXXpu77747dXV1+dznPpeGhoZ3fxQAAAB0OC2adGn06NEZPXp0kuSiiy7a5DZdunRJpVLZ5Lo33ngjDzzwQKZNm5bDDz88SXLVVVflyCOPzLx58zJq1Kht6R0AAIAOrLDvsD711FMZOXJkjjvuuHz5y1/O8uXLm9YtWLAg69atyxFHHNG0bLfddss+++yT+fPnF9UCAAAAHUghj7UZNWpUjjnmmOy5555ZvHhxrrvuupx++ul58MEH06VLlyxbtizV1dXp3bt3s/f17ds3y5Yte8f61dVVqa3tXkSr/7/eToXWg/ZiLJdfQ/72OBHoiJx/mnNOpqMwlimTQgLr2LFjm/49ePDg7LfffhkzZkx++ctf5thjj93s+xobG1tUv6GhMStX1r/rPt9SW9u90HrQXozl8qtUehb67Mvkb8+/hDJw/mnOOZmOwlimPVQqPTe5vFUea9O/f//0798/ixYtSpL069cvDQ0NefXVV5ttt2LFivTr1681WgAAAGA71yqBdcWKFVmyZEl23XXXJMn++++fzp0757HHHmva5pVXXskLL7yQYcOGtUYLAAAAbOdadEtwXV1dXnzxxSTJhg0b8vLLL+fZZ5/NLrvskl122SXXX399jj322FQqlSxevDjXXntt+vTpk6OPPjpJ0rNnz5x00km56qqr0rdv39TW1mb69OkZPHhwDjvssNY7OgAAALZbLQqsCxYsyPjx45tez5o1K7Nmzcq4cePy1a9+Nc8991zmzJmTN954I5VKJSNGjMh1112XHj16NL3n4osvTqdOnXL++ednzZo1GTlyZK666qpUV5uMBAAAgI21KLCOGDEiCxcu3Oz622+//R1rdO3aNZdeemkuvfTSlncHAADADqtVvsMKAAAA75bACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEAptSiwPvHEEzn33HMzatSoDB48OA8++GDTunXr1uXqq6/Oxz72sRx44IE54ogjcsEFF+Tll19uVuO0007L4MGDm/2cf/75xR4NAAAAHUanlmxUX1+fQYMG5ROf+EQuvPDCZuvWrFmTZ555Jp///Ofzvve9L6tWrcqMGTNy1lln5Uc/+lE6dfr7Lk488cRMnDix6XVNTU1BhwEAAEBH06LAOnr06IwePTpJctFFFzVb17Nnz9x5553Nll1++eUZO3ZsXnjhhQwePLhpebdu3VKpVN5tzwAAAOwAWuU7rKtWrUqS7LLLLs2Wz507NyNGjMjYsWMzc+bMpu0AAADg/2rRFdat8eabb2bGjBk58sgj8573vKdp+QknnJDdd989u+66a55//vl8/etfzx/+8IeNrs5uSnV1VWpruxfWY3X1ToXWg/ZiLAPtyfmnOedkOgpjmTIpNLCuX78+kyZNyhtvvJHZs2c3W/fP//zPTf8ePHhw9tprr5x88sn5/e9/n/3222+LdRsaGrNyZX1hfdbWdi+0HrQXY7n8KpWe7d0CtBrnn+ack+kojGXaw+Y+MxV2S/D69eszceLELFy4MN/61rfSu3fvLW6///77p7q6On/+85+LagEAAIAOpJArrOvWrcvEiRPz3HPP5a677mrRxErPPfdcGhoaTMIEAADAJrUosNbV1eXFF19MkmzYsCEvv/xynn322eyyyy7Zdddd88UvfjFPP/10brrpplRVVWXp0qVJ/jaDcE1NTV588cX86Ec/yujRo9O7d++88MILmTFjRj7wgQ9k+PDhrXd0AAAAbLdaFFgXLFiQ8ePHN72eNWtWZs2alXHjxuW8887Lz372syR/e87q202fPj0nnnhiOnfunP/8z//MXXfdlbq6uuy2224ZPXp0zjvvvFRXVxd4OAAAAHQULQqsI0aMyMKFCze7fkvrkmS33XbLd7/73a3rDAAAgB1a4Y+1AQA6vjXrGgqdBXv12vVZ9frqwuoB0DEIrADAVqvpXJ0BU+YWVm/RjLFZVVg1ADqKwh5rAwAAAEUSWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUurU3g0AlEmPXt3SratTIwBAGfhUBvA23bp2yoApcwurt2jG2MJqAQDsaNwSDAAAQCkJrAAAAJSSwAoAAEAp+Q4rANDu1qxrSKXSs9Caq9euz6rXVxdaE4C2JbACAO2upnN1oROeJX+b9GxVoRUBaGtuCQYAAKCUBFYAAABKqUWB9Yknnsi5556bUaNGZfDgwXnwwQebrW9sbMysWbNyxBFHZMiQITnttNPyxz/+sdk2r732WiZNmpSDDjooBx10UCZNmpTXX3+9uCMBAACgQ2lRYK2vr8+gQYNyySWXpKamZqP1t956a+64445ceumluf/++9OnT5+cccYZWbXq798cueCCC/LMM8/k1ltvzW233ZZnnnkmkydPLu5IAAAA6FBaFFhHjx6diRMn5iMf+Uh22qn5WxobG/Od73wn55xzTo477rgMGjQoM2fOTF1dXR5++OEkyQsvvJBf/epXufzyyzN8+PAMGzYsX/va1/KLX/wi//M//1P8UQEAALDde9ezBP/lL3/J0qVLc/jhhzctq6mpySGHHJL58+fn1FNPzfz589O9e/cMHz68aZuDDjoo3bt3z/z58zNw4MAt7qO6uiq1td3fbatvq7dTofWgvRjLAFvWludI52Q6CmOZMnnXgXXp0qVJkn79+jVb3rdv3yxZsiRJsmzZsvTp0ydVVVVN66uqqtKnT58sW7bsHffR0NCYlSvr322rTWpruxdaD9qLsVy8op8DCbSvtjxHOifTURjLtIfNfQYrbJbgt4fRlq5vbGx8x/cBAACwY3rXgbVSqST5+5XWtyxfvrzpqmu/fv2yfPnyNDY2Nq1vbGzMq6++mr59+77bFgAAAOiA3nVg3XPPPVOpVDJv3rymZWvXrs2TTz6ZYcOGJUmGDRuW+vr6zJ8/v2mb+fPnp76+vmkbAAAAeLsWfYe1rq4uL774YpJkw4YNefnll/Pss89ml112ye67757x48fnpptuysCBAzNgwIDMnj073bt3zwknnJAk2WeffTJq1KhcdtllueKKK9LY2JjLLrssRx555DtOuAQAAMCOqUWBdcGCBRk/fnzT61mzZmXWrFkZN25cZsyYkbPPPjtr167N5Zdfntdeey1Dhw7NHXfckR49ejS955prrsnUqVNz5plnJknGjBmTr3zlKwUfDgAAAB1FiwLriBEjsnDhws2ur6qqyoQJEzJhwoTNblNbW5trrrlm6zsEAABgh1TYLMEAAABQJIEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUurU3g0AALSGNesaUqn0LKze6rXrs+r11YXVA+CdCawAQIdU07k6A6bMLazeohljs6qwagC0hFuCAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCXPYQW2az16dUu3rk5lAAAdkU95wHatW9dOGTBlbmH1Fs0YW1gtAADeHbcEAwAAUEoCKwAAAKUksAIAAFBKAisAAAClJLACAABQSgIrAAAApSSwAgAAUEoCKwAAAKXUqYgiY8aMyeLFizdaPnr06Nxyyy2ZNWtWrr/++mbr+vXrl8cee6yI3QMAANABFRJY77///jQ0NDS9Xrp0aU488cR89KMfbVq2995756677mp6XV1dXcSuAQAA6KAKCax9+vRp9vr+++9Pjx498pGPfOTvO+rUKZVKpYjdAQAAsAMoJLC+XWNjY+6///780z/9U7p169a0/KWXXsqoUaPSuXPnDB06NBMnTsxee+1V9O4BAADoIAoPrI899lj+8pe/5OSTT25aNmTIkEyfPj0DBw7MihUrMnv27Jx66ql5+OGH07t373esWV1dldra7oX1WF29U6H1oL0YywBta0vnXOdkOgpjmTIpPLDed999OeCAA/L+97+/adno0aObbTN06NAcffTRmTNnTs4444x3rNnQ0JiVK+sL67G2tnuh9aC9GMtJpdKzvVsAdiBbOuc6J9NRGMu0h819piv0sTbLly/Pz3/+85xyyilb3G7nnXfOvvvum0WLFhW5ewAAADqQQgPrgw8+mM6dO+f444/f4nZr167Nn/70J5MwAQAAsFmF3RL81mRLY8eOTY8ePZqtmzlzZo488sjstttuWbFiRW688cbU19dn3LhxRe0eAACADqawwPpf//VfWbRoUa6++uqN1r3yyiuZOHFiVq5cmd69e+fAAw/Mfffdlz322KOo3QMAtKo16xre8XvzW/O9+tVr12fV66vfbVsAHVphgfXQQw/NwoULN7nuG9/4RlG7AQBoFzWdqzNgytzC6i2aMTarCqsG0DEV+h1WAAAAKIrACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApdWrvBoAdR49e3dKtq9MOAAAt45Mj0Ga6de2UAVPmFlpz0YyxhdYDAKA83BIMAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUksAKAABAKQmsAAAAlFKn9m4AAGBHtGZdQyqVnoXVW712fVa9vrqwegBlILACALSDms7VGTBlbmH1Fs0Ym1WFVQMoh0JuCZ41a1YGDx7c7Ofwww9vWt/Y2JhZs2bliCOOyJAhQ3Laaaflj3/8YxG7BgAAoIMq7Arr3nvvnbvuuqvpdXV1ddO/b7311txxxx2ZMWNG9t5779xwww0544wz8uMf/zg9evQoqgUAAAA6kMImXerUqVMqlUrTT58+fZL87erqd77znZxzzjk57rjjMmjQoMycOTN1dXV5+OGHi9o9AAAAHUxhgfWll17KqFGjMmbMmJx//vl56aWXkiR/+ctfsnTp0ma3CNfU1OSQQw7J/Pnzi9o9AAAAHUwhtwQPGTIk06dPz8CBA7NixYrMnj07p556ah5++OEsXbo0SdKvX79m7+nbt2+WLFnSovrV1VWpre1eRKv/v95OhdaD9mIsA/B2fidQBJ8vKJNCAuvo0aObvR46dGiOPvrozJkzJ0OHDk2SVFVVbXP9hobGrFxZ/656fLva2u6F1oP2sr2N5SIf3wDAxran3wmU1/b2+YKOYXOfEwu7Jfjtdt555+y7775ZtGhRKpVKkjRdaX3L8uXLN7rqCgAAAG9plcC6du3a/OlPf0qlUsmee+6ZSqWSefPmNVv/5JNPZtiwYa2xewAAADqAQm4JnjlzZo488sjstttuWbFiRW688cbU19dn3Lhxqaqqyvjx43PTTTdl4MCBGTBgQGbPnp3u3bvnhBNOKGL3AAAAdECFBNZXXnklEydOzMqVK9O7d+8ceOCBue+++7LHHnskSc4+++ysXbs2l19+eV577bUMHTo0d9xxh2ewAgAAsFmFBNZvfOMbW1xfVVWVCRMmZMKECUXsDgAAgB1Aq3yHFQAAAN4tgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUurU3g0AAPDurVnXkEqlZ6E1V69dn1Wvry60JsDWEFgBADqAms7VGTBlbqE1F80Ym1WFVgTYOm4JBgAAoJQEVgAAAEpJYAUAAKCUBFYAAABKSWAFAACglARWAAAASklgBQAAoJQEVgAAAEqpU3s3ABSnR69u6da1uP/Wq9euz6rXVxdWDwAAtobACh1It66dMmDK3MLqLZoxNqsKqwYAAFvHLcEAAACUksAKAABAKbklGACANmGuBWBrCawAALQJcy0AW0tgBTZrzbqGVCo927sNAAB2UAIrsFk1nasL/0s4AAC0VCGB9eabb85PfvKT/OlPf0qXLl1y4IEHZuLEiRk0aFDTNlOmTMkPfvCDZu8bOnRo7rvvviJaAAAAoIMpJLA+/vjj+dSnPpUDDjggjY2N+dd//decccYZmTt3bmpra5u2O+ywwzyAud0AAAw3SURBVHLVVVc1ve7cuXMRuwcAAKADKiSw3n777c1eX3XVVTn44IPz29/+NmPGjGla3qVLl1QqlSJ2CQAAQAfXKt9hraury4YNG9KrV69my5966qmMHDkyvXr1yiGHHJLzzz8/ffv2bY0WAAAA2M61SmC98sor8/73vz/Dhg1rWjZq1Kgcc8wx2XPPPbN48eJcd911Of300/Pggw+mS5cuW6xXXV2V2truhfVXXb1TofWgvRjLALSm7WG2eL8Hi+fzBWVSeGCdPn16nnrqqdx7772prq5uWj527N9nBx08eHD222+/jBkzJr/85S9z7LHHbrFmQ0NjVq6sL6zH2truhdaD9vJ/x3LZP1QAsH3ZHmaL95mueD4r0x429zm20MA6bdq0PPLII/n2t7+dvfbaa4vb9u/fP/3798+iRYuKbAEAAIAOorDAOnXq1DzyyCO56667ss8++7zj9itWrMiSJUuy6667FtUCAAAAHUghgfVrX/tafvjDH+aGG25Ir169snTp0iRJ9+7ds/POO6euri7XX399jj322FQqlSxevDjXXntt+vTpk6OPPrqIFgAAAOhgCgms99xzT5LkM5/5TLPl5513XiZMmJDq6uo899xzmTNnTt54441UKpWMGDEi1113XXr06FFECwAAAHQwhQTWhQsXbnF9TU3NRs9qBQCAd6M1ZjFevXZ9Vr2+utCawLZrlcfaAABAayt6FuPkbzMZryq0IvBu7NTeDQAAAMCmCKwAAACUksAKAABAKQmsAAAAlJLACgAAQCkJrAAAAJSSwAoAAEApCawAAACUUqf2bgB2ZD16dUu3ru/uv2Gl0rOgbgCANesaCv3dumZdQ2o6VxdWb/Xa9Vn1+urC6kHZCazQjrp17ZQBU+YWVm/RjLGF1QKAHVFN5+rCfzcXXW9VYdWg/NwSDAAAQCkJrAAAAJSSwAoAAEApCawAAACUkkmXoIWKmNEXAABoOZ++oYWKntE3MasvAABsiVuCAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJbME02F5DA0A0NGsWdeQSqVnYfVWr12fVa+vLqweFM2neTqsoh9D4xE0AEB7q+lcXfjnm1WFVYPiCawAAEBhir7LzVXgHZvACgAAFKY17nJzFXjHZdIlAAAASskVVgAA2EFtbhKnIid2gnejzQPr3Xffndtvvz1Lly7NP/7jP+biiy/OwQcf3NZtAADADq/oSZwSE1VSrDYNrI888kimTZuWyy67LAcddFDuueeenH322Zk7d2523333tmwFAACgECaaaj1tGljvvPPOjBs3LqecckqS5NJLL82vfvWr3HvvvbngggvaspXC7WiDtDWecbpmXUNqOlcXWhMAgO3b9vDs2aInmvrDFR8p/TG3lTYLrG+++WZ+//vf58wzz2y2/PDDD8/8+fPbqo1Ws6PNhlb08SZ/O2bPTQUA4O12xGfP7ojHvDlVjY2NjW2xo7/+9a/50Ic+lO9+97s55JBDmpZff/31eeihh/Loo4+2RRsAAABsJ9r8sTZVVVUtWgYAAMCOrc0Ca+/evVNdXZ2lS5c2W758+fL069evrdoAAABgO9FmgbVLly7Zb7/9Mm/evGbL582bl2HDhrVVGwAAAGwn2nSW4DPOOCOTJ0/OkCFDMnz48Nx7771ZsmRJTj311LZsAwAAgO1AmwbW448/Pq+++mpmz56dJUuWZNCgQbnllluyxx57tGUbAAAAbAfabJZgAAAA2BptPkswAAAAtESHDKx33313xowZkwMOOCAnnnhinnzyyS1u//jjj+fEE0/MAQcckKOOOir33ntvG3UKm7c14/gnP/lJzjzzzBx66KEZNmxYTj755PzsZz9rw25h87b2nPyWJ598Mh/4wAdywgkntHKH0DJbO5bffPPNfPOb38yYMWOy//7758Mf/nC+853vtFG3sHlbO5YfeuihfPzjH8/QoUNz+OGH50tf+tJGT/6A1tLhAusjjzySadOm5dxzz82cOXMybNiwnH322Xn55Zc3uf1LL72Uc845J8OGDcucOXPyuc99LlOnTs2jjz7axp3D323tOH788cdz6KGH5pZbbsmcOXMyevTonHfeeS0OBtBatnYsv+W1117LhRdemJEjR7ZRp7Bl2zKWL7jggvzqV7/KFVdckR//+Mf55je/mcGDB7dh17CxrR3LTz31VCZPnpxx48bl4Ycfzg033JAXXnghX/rSl9q4c3ZUHe47rCeffHIGDx6cqVOnNi079thjc9xxx+WCCy7YaPurr746P/3pT/OTn/ykadkll1yS559/Pt///vfbpGf4v7Z2HG/KJz/5yRx88MGZMmVKa7UJ72hbx/J5552X973vfWlsbMyjjz6ahx9+uC3ahc3a2rH861//Ol/84hfz05/+NH369GnLVmGLtnYs33777fnud7+bX/ziF03LHnjggUydOjXz589vk57ZsXWoK6xvvvlmfv/73+fwww9vtvzwww/f7H+o3/3udxttf8QRR2TBggVZt25dq/UKm7Mt43hT6urq0qtXr6Lbgxbb1rF89913Z9myZfn85z/f2i1Ci2zLWP73f//3HHDAAfnWt76VD33oQzn22GMzderU1NXVtUXLsEnbMpaHDx+epUuX5uc//3kaGxuzYsWKPPLII/nQhz7UFi1Dxwqsr776ahoaGtKvX79my/v27bvZ++yXLVuWvn37NlvWr1+/rF+/Pq+++mqr9Qqbsy3j+P+6++6788orr+TjH/94a7QILbItY3nhwoW54YYbcvXVV6e6urot2oR3tC1j+aWXXspTTz2VP/zhD5k1a1YuvfTS/OpXv8pFF13UFi3DJm3LWB42bFi+/vWv50tf+lL233//jBw5Mo2NjZk5c2ZbtAwdK7C+paqqqkXLNrfurbukt/QeaG1bO47f8uijj+aqq67KNddc4xnHlEJLx/Kbb76ZiRMnZvLkydlrr73aojXYKltzXm5sbExVVVW+/vWvZ+jQoRk1alQuvfTSPProo1m2bFlrtwpbtDVj+fnnn8/UqVPzL//yL3nggQdy2223ZenSpfnKV77S2m1CkqRTezdQpN69e6e6unqjvxAtX758o78kvaVfv34b/eJYvnx5OnXqlNra2lbrFTZnW8bxWx599NFMnjw5M2fOzFFHHdWabcI72tqxvGTJkjz//PO5+OKLc/HFFydJNmzYkMbGxnzgAx/ILbfckiOOOKJNeoe325bzcqVSSf/+/dOzZ8+mZfvss0+S5OWXX37H8zm0hm0ZyzfffHOGDBmSs846K0nyvve9L926dcunP/3pnH/++dltt91avW92bB3qCmuXLl2y3377Zd68ec2Wz5s3L8OGDdvkew488MBNbr///vunc+fOrdYrbM62jOPkb7P+TZo0KdOnT89HPvKR1m4T3tHWjuX+/fvnoYceypw5c5p+Tj311Lz3ve9tmskS2sO2nJeHDx+eJUuWNPvO6qJFi5LE3S+0m20Zy2vWrNnoKxpvve5gc7dSUtVf/epXv9reTRSpR48emTVrViqVSmpqanLjjTfmySefzLRp09KrV69Mnjw5P/3pT3PMMcckSf7hH/4ht956a5YvX5499tgjP/vZz3LTTTdlypQp2Xfffdv5aNhRbe04njt3biZPnpxJkyblwx/+cOrr61NfX59169alpqamnY+GHdnWjOXq6ur07du32c/TTz+dRYsWZcKECenSpUt7Hw47sK09L++999554IEHsmDBguy7775ZtGhRLr/88hxyyCE56aST2vlo2JFt7Vhes2ZNbr/99vTu3Tu1tbV5/vnnc+WVV6ZSqTRddYXW1KFuCU6S448/Pq+++mpmz56dJUuWZNCgQbnlllua/pr5v//7v82232uvvXLLLbdk+vTpuffee7PrrrvmkksuyXHHHdce7UOSrR/H3/ve97J+/fpMmzYt06ZNa1r+wQ9+MHfddVeb9g5vt7VjGcpqa8fyzjvvnDvvvDNTp07NJz/5yfTq1StHH310ix9NBq1la8fyiSeemLq6utx9992ZOXNmevbsmREjRmTSpEnt0T47oA73HFYAAAA6hg71HVYAAAA6DoEVAACAUhJYAQAAKCWBFQAAgFISWAEAACglgRUAAIBSElgBAAAoJYEVAACAUhJYAQAAKKX/B5n1RgcvBBlKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Reinitialize using the same entropy\n",
    "rs = np.random.RandomState(entropy)\n",
    "\n",
    "bs = StationaryBootstrap(opt.loc[\"Mkt-RF\", \"stationary\"], excess_market, random_state=rs)\n",
    "results = bs.apply(sharpe_ratio, 2500)\n",
    "SR = pd.DataFrame(results[:, -1:], columns=['SR'])\n",
    "fig = SR.hist(bins=40)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probit (statsmodels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second example makes use of a Probit model from statsmodels.  The demo data is university admissions data which contains a binary variable for being admitted, GRE score, GPA score and quartile rank. This data is downloaded from the internet and imported using pandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            admit         gre         gpa       rank\n",
      "count  400.000000  400.000000  400.000000  400.00000\n",
      "mean     0.317500  587.700000    3.389900    2.48500\n",
      "std      0.466087  115.516536    0.380567    0.94446\n",
      "min      0.000000  220.000000    2.260000    1.00000\n",
      "25%      0.000000  520.000000    3.130000    2.00000\n",
      "50%      0.000000  580.000000    3.395000    2.00000\n",
      "75%      1.000000  660.000000    3.670000    3.00000\n",
      "max      1.000000  800.000000    4.000000    4.00000\n"
     ]
    }
   ],
   "source": [
    "import arch.data.binary\n",
    "\n",
    "binary = arch.data.binary.load()\n",
    "binary = binary.dropna()\n",
    "print(binary.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting the model directly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first steps are to build the regressor and the dependent variable arrays.  Then, using these arrays, the model can be estimated by calling `fit`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Const   -3.003536\n",
      "gre      0.001643\n",
      "gpa      0.454575\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "\n",
    "endog = binary[['admit']]\n",
    "exog = binary[['gre', 'gpa']]\n",
    "const = pd.Series(np.ones(exog.shape[0]), index=endog.index)\n",
    "const.name = 'Const'\n",
    "exog = pd.DataFrame([const, exog.gre, exog.gpa]).T\n",
    "\n",
    "# Estimate the model\n",
    "mod = sm.Probit(endog, exog)\n",
    "fit = mod.fit(disp=0)\n",
    "params = fit.params\n",
    "print(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The wrapper function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most models in statsmodels are implemented as classes, require an explicit call to `fit` and return a class containing parameter estimates and other quantities.  These classes cannot be directly used with the bootstrap methods.  However, a simple wrapper can be written that takes the data as the only inputs and returns parameters estimated using a Statsmodel model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def probit_wrap(endog, exog):\n",
    "    return sm.Probit(endog, exog).fit(disp=0).params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A call to this function should return the same parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Const   -3.003536\n",
       "gre      0.001643\n",
       "gpa      0.454575\n",
       "dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "probit_wrap(endog, exog)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The wrapper can be directly used to estimate the parameter covariance or to construct confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          Const           gre       gpa\n",
      "Const  0.384854 -5.813301e-05 -0.101077\n",
      "gre   -0.000058  4.244940e-07 -0.000057\n",
      "gpa   -0.101077 -5.701083e-05  0.039506\n"
     ]
    }
   ],
   "source": [
    "from arch.bootstrap import IIDBootstrap\n",
    "\n",
    "bs = IIDBootstrap(endog=endog, exog=exog)\n",
    "cov = bs.cov(probit_wrap, 1000)\n",
    "cov = pd.DataFrame(cov, index=exog.columns, columns=exog.columns)\n",
    "print(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Const    0.620366\n",
      "gre      0.000652\n",
      "gpa      0.198760\n",
      "dtype: float64\n",
      "T-stats\n",
      "Const   -4.841552\n",
      "gre      2.521038\n",
      "gpa      2.287054\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "se = pd.Series(np.sqrt(np.diag(cov)), index=exog.columns)\n",
    "print(se)\n",
    "print('T-stats')\n",
    "print(params / se)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          Const       gre       gpa\n",
      "Lower -4.194772  0.000354  0.020491\n",
      "Upper -1.582297  0.002943  0.825360\n"
     ]
    }
   ],
   "source": [
    "ci = bs.conf_int(probit_wrap, 1000, method='basic')\n",
    "ci = pd.DataFrame(ci, index=['Lower', 'Upper'], columns=exog.columns)\n",
    "print(ci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Speeding things up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting values can be provided to `fit` which can save time finding starting values.  Since the bootstrap parameter estimates should be close to the original sample estimates, the full sample estimated parameters are reasonable starting values.  These can be passed using the `extra_kwargs` dictionary to a modified wrapper that will accept a keyword argument containing starting values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def probit_wrap_start_params(endog, exog, start_params=None):\n",
    "    return sm.Probit(endog, exog).fit(start_params=start_params, disp=0).params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          Const           gre       gpa\n",
      "Const  0.384854 -5.813301e-05 -0.101077\n",
      "gre   -0.000058  4.244940e-07 -0.000057\n",
      "gpa   -0.101077 -5.701083e-05  0.039506\n"
     ]
    }
   ],
   "source": [
    "bs.reset()  # Reset to original state for comparability\n",
    "cov = bs.cov(\n",
    "    probit_wrap_start_params,\n",
    "    1000,\n",
    "    extra_kwargs={'start_params': params.values})\n",
    "cov = pd.DataFrame(cov, index=exog.columns, columns=exog.columns)\n",
    "print(cov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrapping Uneven Length Samples\n",
    "Independent samples of uneven length are common in experiment settings, e.g., A/B testing of a website.  The `IIDBootstrap` allows for arbitrary dependence within an observation index and so cannot be naturally applied to these data sets. The `IndependentSamplesBootstrap` allows datasets with variables of different lengths to be sampled by exploiting the independence of the values to separately bootstrap each component. Below is an example showing how a confidence interval can be constructed for the difference in means of two groups. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from arch.bootstrap import IndependentSamplesBootstrap\n",
    "\n",
    "\n",
    "def mean_diff(x, y):\n",
    "    return x.mean() - y.mean()\n",
    "\n",
    "\n",
    "rs = np.random.RandomState(0)\n",
    "treatment = 0.2 + rs.standard_normal(200)\n",
    "control = rs.standard_normal(800)\n",
    "\n",
    "bs = IndependentSamplesBootstrap(treatment, control, random_state=rs)\n",
    "print(bs.conf_int(mean_diff, method='studentized'))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
